/* Copyright (c) CERN
 * * Artistic License 2.0
 * */

import java.util.Random;
import java.lang.Math;
import org.apache.commons.math.special.Erf;
import java.util.Arrays;
import org.apache.commons.math.MathException;

/**
 * HalfNormal distribution.
 * This is not the classic half-normal distribution, i.e. the absolute value of a
 * centered normal random variable. Instead, this is the right part of a normal random variable
 * of standard deviation 'sigma' and centered at 'max'.
 *
 * @author Xavier Grehant
 */
public class HalfNormal extends DistroUtils implements Distro {

/**
 * Distribution parameter.
 * All values are lower than 'max', distributed as if they were the values of
 * a gaussian centered around 'max'.
 */
   int max;

/**
 * Distribution parameter.
 * The standard deviation of a gaussian centered at 'max', whose part at the right
 * of 'max' fits this distribution.
 */
   double sigma;

/**
 * A uniform number generator.
 */
   Random rand;

/**
 * Constructor taking distribution parameters.
 * @param max center of the corresponding gaussian.
 * @param sigma standard deviation of the corresponding gaussian.
 */
   public HalfNormal(int max, double sigma) {
      construct(max, sigma);
   }

/**
 * Constructor taking a set of samples.
 * @param array array of samples.
 */
   public HalfNormal(int[] array) {
      construct(array);
   }

/**
 * This method transforms samples into samples of the corresponding gaussian.
 * @return an array of integers following the corresponding gaussian, of twice
 *         the input size.
 * @param array sorted samples of the distribution.
 */
   int[] mirror(int[] array) {
      // shift to the left of zero axis and mirror.
      max = array[array.length - 1];
      int[] left = addAll(array, - max);
      int[] right = multiplyAll(left, - 1);
      left = flip(left); // just for viewing the bell curve
      return concat(left, Arrays.copyOfRange(right, 1, right.length - 1));
   }

/**
 * This method builds a sample set including samples of the input sets.
 * @return set of integers containing integers of the input sets.
 * @param left set of integers.
 * @param right set of integers.
 */
   int[] concat(int[] left, int[] right) {
      int[] result = new int[left.length + right.length];
      for (int i = 0; i < left.length; i++) result[i] = left[i];
      for (int i = 0; i < right.length; i++) result[left.length + i] = right[i];
      return result;
   }

/**
 * Adds the same number to each number of an array.
 * @return the array after addition on all its members.
 * @param array the array
 * @param addend the addend
 */
   int[] addAll(int[] array, int addend) {
      for (int i = 0; i < array.length; i++) array[i] += addend;
      return array;
   }
   
/**
 * Multiplies by the same multiplier each number of an array.
 * @return the array after multiplication on all its members.
 * @param array the array.
 * @param multiplier the multiplier.
 */
   int[] multiplyAll(int[] array, int multiplier) {
      for (int i = 0; i < array.length; i++) array[i] *= multiplier;
      return array;
   }

/**
 * Constructor that takes a set of samples and discards lower values
 * for the estimation.
 * @param array the set of samples
 * @param minvalue the minimum value of interest
 */
   public HalfNormal(int[] array, int minvalue) {
      int size = 0;
      for (int i = 0; i < array.length; i++)
         if (array[i] >= minvalue) size++;
      if (size == 0) construct(0, 0f);
      else {
         int[] result = new int[size];
         int index = 0;
         for (int i = 0; i < array.length; i++) {
            if (array[i] >= minvalue) {
               result[index] = array[i];
               index++;
            }
         }
         construct(result);
      }
   }

/**
 * This method sets instance fields.
 * @return 
 * @param max value to assign to instance field max.
 * @param sigma value to assign to instance field sigma.
 */
   void construct(int max, double sigma) {
      this.max = max;
      this.sigma = sigma;
      rand = new Random();
      System.out.println("max = " + max + ", sigma = " + sigma);
   }

/**
 * This method estimates the fit after discarding some upper samples
 * supposed to be noise. Noise in upper samples affects the distribution
 * badly since the maximum value is a parameter of the distribution.
 * @return 
 * @param array the set of samples.
 */
   void construct(int[] array) {
      array = discard(array, .01, 1);
      // remove less than 1o/oo samples if they have a std of more than 2%
      Arrays.sort(array);
      construct(array[array.length - 1], std(mirror(array), 0d));
   }

/**
 * This factory method estimates a biased GPD to correctly fit samples over minvalue.
 * @return a biased half-normal distribution fitting the samples over 'minvalue'.
 * @param values the sample set.
 * @param minvalue the lowest value of interest.
 * @param refinements the number of times the estimation is biased.
 */
   static HalfNormal biased(int[] values, int minvalue, int refinements) {
      HalfNormal hn = new HalfNormal(values);
      for (int refined = 0; refined < refinements; refined++) {
         values = hn.bias(values, minvalue);
         hn = new HalfNormal(values);
      }
      return hn;
   }

/**
 * Random number generator according to the distribution.
 * @return a random double.
 */
   public double nextDouble() {
      return max - Math.abs(rand.nextGaussian() * sigma);
   }

/**
 * Cumulative distribution function of the distribution.
 * @return P(X <= 'limit')
 * @param limit higher bound.
 */
   public double cdf(int limit) throws MathException {
      if (limit >= max) return 1; else
      return 1d + Erf.erf((double) (limit - max) / ((double) sigma * Math.sqrt(2)));
   }
}
